%% TCM модуляция

% Инициализация
clear all; % удаление всех переменных, функций и ссылок
clc;       % очистка командного окна

numFig = 0; % счетчик графиков

type = 'PSK';

% Кодирование
switch type
    case 'QAM'
        % Описание кода для 8-QAM
        trel = struct('numInputSymbols', 4, 'numOutputSymbols', 8,...
                      'numStates',       4, ...
                      'nextStates', [0 0 1 1; 2 2 3 3; 0 0 1 1; 2 2 3 3],...
                      'outputs',    [0 5 1 4; 3 6 2 7; 1 4 0 5; 2 7 3 6]);
    case 'PSK'
        % Описание кода для 8-PSK
        trel = struct('numInputSymbols', 4, 'numOutputSymbols', 8,...
                      'numStates',       4, ...
                      'nextStates', [0 0 1 1; 2 2 3 3; 0 0 1 1; 2 2 3 3],...
                      'outputs',    [0 4 2 6; 1 5 3 7; 2 6 0 4; 3 7 1 5]);
    otherwise
        disp('error: неизвестный тип модуляции. Известны "QAM" и "PSK"');
end

% Декодирование
K     = log2(trel.numStates) + 1; % кодовое ограничение
tblen = 5 * K;                    % глубина просмотра при декодировании

% Модуляция
q = trel.numOutputSymbols;
m = log2(q);
n = m;
k = log2(trel.numInputSymbols); % число бит информации
R = k / n;

% без кодирования
qUncoding = trel.numInputSymbols;
mUncoding = log2(qUncoding);

E = 1; % энергия сигнала

signals = modulationSignals(type, q, E); % сигнальные точки

% сигнальное множество
numFig = numFig+1;
figure(numFig);
plot(signals(1, :), signals(2, :), 'ko');
for i = 1 : q
    text( signals(1, i), signals(2, i), [' ', num2str(i-1)], 'Color', 'b' );
end
title(['Signal set of coded systems, TCM(k=', num2str(k), ', n=',...
       num2str(n), ', K=', num2str(K), ', ', num2str(q), '-', type, ')'])
xlabel('s_{i1}')
ylabel('s_{i2}')
axis('square');

% без кодирования
[signalsUncoding, dUncoding] = modulationSignals(type, qUncoding,...
                                                 E); % сигнальные точки

% сигнальное множество
numFig = numFig+1;
figure(numFig);
plot(signalsUncoding(1, :), signalsUncoding(2, :), 'k*');

for i = 1 : qUncoding
    text( signalsUncoding(1, i), signalsUncoding(2, i), [' ', num2str(i-1)], 'Color', 'b' );
end
title(['Signal set of uncoded systems, ', num2str(qUncoding), '-', type])
xlabel('s_{i1}')
ylabel('s_{i2}')
axis('square');


%% TCM

% SNRdBinit = 4;   % начальное значение отношения сигнал/шум
% SNRdBfin  = 12;   % конечное значение отношения сигнал/шум
% dSNRdB    = 2;   % шаг отношения сигнал/шум

nErrMax = 20; % число ошибок, при котором завершается моделирование при
               % данном значении отношения сигнал/шум. Происходит переход
               % к следующему значению отношения сигнал/шум
               
listSNRdB = [4, 6, 8]; % , 10, 10.25, 10.5 9, 9.5, 10

% практические вероятности ошибок при различном соотношении сигнал/шум
Pe_bit         = zeros( size(listSNRdB) ); % zeros( size(SNRdBinit:dSNRdB:SNRdBfin) );
Pe             = Pe_bit;
PeUncoding     = Pe_bit;
Pe_bitUncoding = Pe_bit; % без кодирования


j = 1; % счетчик отношений сигнал к шуму
% for SNRdB = SNRdBinit : dSNRdB : SNRdBfin
for SNRdB = listSNRdB
    
    % SNR = E/No
    % SNR - signal to noise ratio (отношение сигнал/шум)
    SNR = 10^(SNRdB/10); % перевод из децибел
    
    noise  = E/SNR;        % шум
    DNoise = noise/2;      % дисперсия шума, для БГШ D = NO/2
    sigma  = sqrt(DNoise); % среднее квадратичное отклонение
    
    l = 0; % счетчик испытаний
    
    % счетчики ошибок
    nBitErr         = 0;
    nErr            = 0;
    nBitErrUncoding = 0;
    nErrUncoding    = 0;
    
    finalstate = []; % состояние решетчатого кодера
    % состояния декодера Витерби
    f1 = []; f2 = []; f3 = [];
    
    msgs = zeros(k*tblen, 1); % история сообщений
    
    while nBitErrUncoding < nErrMax || nBitErr < nErrMax
        % Случайные данные
        msg = randi([0 1], [k, 1]);
        
        % Модуляция TCM
        [point, finalstate] = tcmEncoder(msg, trel, signals, finalstate);
        
        % без кодирования
        msgNum        = bi2de(msg', 'left-msb');
        pointUncoding = signalsUncoding(:, msgNum + 1);
        
        % Добавление шума
        noisePoint    = [sigma*randn; sigma*randn];
        point         = point         + noisePoint;
        pointUncoding = pointUncoding + noisePoint;
        
        % Демодуляция TCM
        [decoded, f1, f2, f3] = tcmDecoder(point, trel, signals, tblen,...
                                           f1, f2, f3);
        
        if l >= tblen && ~isequal(msgs(1:k), decoded)
            % сосчитать число битовых ошибок
            nBitErr = nBitErr + biterr(msgs(1:k), decoded);
            nErr = nErr + 1;
        end
        
        % обеспечить задержку декодирования
        msgs(1:k) = msg;
        msgs      = circshift(msgs, -k); % циклически сдвинуть 
                                         % элементы массива
        
        % без кодирования
        [~, z] = min( pdist2(signalsUncoding', pointUncoding',...
                      'euclidean') ); % Евклидова метрика
        z = z - 1; % переход от позиции сигнала к его номеру
        
        if ~isequal(msgNum, z)
            % сосчитать число битовых ошибок
            nBitErrUncoding = nBitErrUncoding + biterr(msgNum, z);
            nErrUncoding = nErrUncoding + 1;
            fprintf('SNRdB = %d/%d, l = %d, nBitErr = %d/%d, nBitErrUncoding = %d/%d\n',...
                    SNRdB, listSNRdB(end), l, nBitErr, nErrMax, nBitErrUncoding, nErrMax);
        end
        
        l = l+1; % перейти к следующему испытанию
    end
    
    % подсчитать вероятность ошибки в данном испытании
    Pe_bit(j)         = nBitErr / ( (l - tblen) * k );
    Pe(j)             = nErr / (l - tblen);
    Pe_bitUncoding(j) = nBitErrUncoding / (l * k);
    PeUncoding(j)     = nErrUncoding / l;
    
    fprintf('\nSNRdB = %d/%d, l = %d, nBitErr = %d/%d, Pe_bit(j) = %d\n\n',...
            SNRdB, listSNRdB(end), l, nBitErr, nErrMax, Pe_bit(j));

    fprintf('\nSNRdB = %d/%d, l = %d, nBitErrUncoding = %d/%d, Pe_bitUncoding(j) = %d\n\n',...
            SNRdB, listSNRdB(end), l, nBitErrUncoding, nErrMax, Pe_bitUncoding(j));
    
    j = j+1; % перейти к следующему соотношению отношения сигнала к шуму
end

%% Вывод графика зависимости вероятности ошибки от отношения сигнал/шум

df = freeDistanceOfTrellis(trel, signals); % свободное расстояние
                                           % решетчатого кода
G  = 20 * log10(df / dUncoding);

% (E/N0)_bit = (E/N0) / (R * m);
SNRdB_bit         = listSNRdB / (R * m);
% (E/N0)_bit = (E/N0) / m;
SNRdB_bitUncoding = listSNRdB / mUncoding;


listSNRdB_Th = 4 : 0.1 : 20;
listSNRdB_bit_Th         = listSNRdB_Th / (R * m);
listSNRdB_bitUncoding_Th = listSNRdB_Th / mUncoding;

%% Теоретические значения
numFig = numFig+1;
figure(numFig);

SNR        = 10.^( listSNRdB / 10 );
listSNR_Th = 10.^( listSNRdB_Th / 10 );

% SNR    = E/noise      => noise  = E/SNR
% DNoise = noise/2      => DNoise = E / (2*SNR)
% sigma  = sqrt(DNoise) => sigma  = sqrt( E / (2*SNR) )
twoSigma = 2 * sqrt(  E ./ ( 2 * listSNR_Th )  );

thoerPe = Q( df ./ twoSigma );

switch type
    case 'QAM'
        if mod(mUncoding, 2) % если прореженное сигнальное множество
            thoerPeUncoding = exp(-(3*SNR/(2*qUncoding-1))); % верхняя граница
            %   Pe_theor = (1/2)*exp(-(3*SNR/(2*q-1))); % кажется, что это тоже верхняя граница
        else         % если цельное сигнальное множество
            thoerPeUncoding = exp(-(3/2)*SNR/(qUncoding-1)); % верхняя граница
            %   Pe_theor = (1/2)*exp(-(3/2)*SNR/(q-1)); % кажется, что это тоже верхняя граница
        end
        % thoerPeUncodingQAM = 1 - (   1 - 2 * Q(    )   );
    case 'PSK'
        % thoerPe_bitUncoding = ( 2 / m ) * Q(  sqrt( 2 * SNRdB_bitUncoding * m ) * sin(pi/q) );
%         thoerPeUncoding = 2 * Q(  sqrt( 2 * SNR ) * sin(pi/qUncoding)  );
        thoerPeUncoding = 2 * Q(  sqrt( 2 * listSNR_Th ) * sin(pi/qUncoding)  );
    otherwise
        disp('error: неизвестный тип модуляции. Известны "QAM" и "PSK"');
end

semilogy(listSNRdB,   PeUncoding,      'k--*',...
         listSNRdB/R, Pe,              'k-o',...
         listSNRdB_Th,   thoerPeUncoding, 'k--',...
         listSNRdB_Th/R, thoerPe,         'k-')
title(['Theoretic P_{e}(E/N_{o})', '; G = ', num2str(G)])
xlabel('E/N_{o}, dB')
ylabel('P_{e}')
legend([num2str(qUncoding), '-', type],...
       ['TCM(', num2str(q), '-', type, ')'], 0)
grid on

thoerPe_bit = thoerPe / k;
thoerPe_bitUncoding = thoerPeUncoding / m;

%% Практические значения
numFig = numFig+1;
figure(numFig);

% график зависимости вероятности ошибки от отношения сигнал/шум
% semilogy(SNRdB_bit, Pe_bit, 'k-o')
semilogy(SNRdB_bitUncoding, Pe_bitUncoding,      'k--*',...
         SNRdB_bit,         Pe_bit,              'k-o',...
         listSNRdB_bitUncoding_Th, thoerPe_bitUncoding, 'k--',...
         listSNRdB_bit_Th,  thoerPe_bit,         'k-')
title(['P_{e}^{(bit)}(E/N_{o})_{(bit)}', '; G = ', num2str(G)])
xlabel('(E/N_{o})_{(bit)}, dB')
ylabel('P_{e}^{(bit)}')
legend([num2str(qUncoding), '-', type],...
       ['TCM(', num2str(q), '-', type, ')'], 0)
grid on